WM-99-118 



O 



Nonperturbative evaluation of the few-body 



O '. states for scalar % interaction 

(N 

Oh' 

Qetin gavkliQ 

^ ! Department of Physics, College of William and Mary, Williamsburg, VA, 

(^ ; 23187 

O 

o 

^ ' Abstract 

' , A knowledge of nonpertubative propagators is often needed when 

Q I the standard perturbative methods are not apphcable. An example of 

1-^ ' this is the bound state problem in field theory. While a nonperturba- 

tive result is valuable by itself, it is also an important guide for those 
K^ \ who work on developing phenomenological models for the nonper- 

;_( ■ turbative problem. The Feynman-Schwinger representation approach 

^ provides a convenient framework for calculating nonperturbative prop- 

agators. In this paper we provide an algorithm for computing 1,2, and 
3 body bound states with the inclusion of all self energies, vertex cor- 
rections, ladder and crossed ladder exchanges. The calculation is done 
in the quenched approximation by ignoring the matter loops. We pro- 
vide simulation results for 1,2 and 3-body states. 

PACS codes: 12.38.-t 

Keywords: Nonperturbative, Monte-Carlo, bound states, Feynman-Schwinger 

representation 



^ csavkli@physics . wm.edu 



Program Summary 

Title of the Programs: phi3 

Computer: Sun 19 

Operating system: Unix 

Programming language used: FORTRAN 77 

Peripherals used: Laser Printer 

Number of lines in distributed program: 1510 

Keywords: Nonperturbative, bound states, Feynman-Schwinger representa- 
tion, Monte-Carlo, numerical evaluation. 

Nature of physical problem: 

The program provided here evaluates the mass and distribution probabilities 
of the fully interacting n-body propagator (n < 3) for scalar x^0 interaction 
in 3+1 dimensions. The evaluation takes into account all self energy, vertex 
dressing, and exchange interaction contributions except those involving mat- 
ter loops (the quenched approximation). 

Method of solution: 

The Feynman-Schwinger representation approach is used to express the field 
theoretical Green's function in terms of a quantum mechanical path integral. 
The resultant expression is evaluated using a Monte-Carlo simulation. 

Restrictions of the program: 

Only n = 1,2, and 3 body propagators are considered. The extension to 

n > 4 requires straightforward modifications in the program. 

Typical running time: 

About 1 day for the 1-body propagator with the self energy. 



LONG WRITE-UP 

1 Introduction 

In nuclear physics one is often faced by problems that require nonperturbative 
methods. The best known example is the problem of bound states. Even if 
the underlying theory may have a small coupling constant (such as in QED), 
and therefore allows the use of perturbation theory in general, the treatment 
of bound states are inherently nonperturbative. The n-body bound state is 
defined by the pole of the interacting n-body propagator. A perturbative 
approximation of n-body propagator does not produce the bound state pole 
location. Therefore it is essential that reliable nonperturbative methods that 
take all orders of interaction into account are developed. For this reason, 
numerous nonperturbative methods have been developed and successfully 
used in the literature. Some of the best known examples are lattice gauge 
theory |I[] (LGT), and relativistic bound state equations 0, |^, ^. 

In a recent paper [^] we (along with authors J. Tjon, and F. Gross) have 
discussed yet another method known as the Feynman-Schwinger representa- 
tion(FSR). The basic idea in the FSR approach 0, H 1^, § is to integrate out 
all fields at the expense of introducing quantum mechanical path integrals 
over the trajectories of particles. Replacing the path integrals over fields with 
path integrals over trajectories has an enormous computational advantage. 
The advantage is due to the fact that the path integration over trajecto- 
ries involves a variation of lines rather than fields in a volume. Therefore 
the degrees of freedom is considerably fewer. The FSR approach differs from 
the LGT in that it utilizes a space-time continuum, therefore maintaining the 
Poincare symmetry. On the other hand it should be pointed out that the FSR 
approach is not without its drawbacks. In particular, how to extend the FSR 
approach to include fermions is not known. In the past researchers P, ^, |10 



have attacked the fermion problem using various approximations. An exact 
result involving fermions is an important problem and requires further study. 
While being able to calculate a nonperturbative result by itself is interest- 
ing, an additional motivation in studying the FSR approach is to determine 
which subsets of diagrams give the dominant contribution to the n-body 
propagator. This is particularly important in determining what kind of ap- 
proximations are reasonable within the context of bound state equations. An 



example of how the FSR results can be used to compare different nonpertur- 
bative approximation schemes was presented in Refs [^, 11|. In those works 



the emphasis was on the development of the formalism and application to the 
1 and 2-body propagators. In Ref |Tl[, The FSR prediction for the 1-body 



mass in Scalar QED was compared by the rainbow-Dyson-Schwinger equa- 
tion prediction. It was found |jll[ that while the FSR approach provides a 



real mass pole for all coupling strengths, the Dyson-Schwinger equation pro- 
vides a complex mass pole beyond a critical coupling strength. Furthermore 
it was found that, for Scalar QED in 0+1 dimension, the vertex corrections 
to the exchange interaction do not contribute to the 2-body binding en- 
ergy 0]. These examples demonstrate the potential usefulness of the FSR 
approach. The knowledge about the nonperturbative propagators and ver- 
tices provided by the FSR is valuable as an input in testing and improving 
the modeling of other nonperturbative approaches such as Dyson-Schwinger 



equations. 12, 13 



The two and three-body bound state sectors provide possibilities for 
the application of the FSR formalism. In particular it is important to see 
how various bound state equations (such as Bethe-Salpeter, Gross (spec- 
tator), Blankenbekler-Sugar, Equal-time, and nonrelativistic) compare with 
the quenched FSR results. Applications of the FSR approach to 2 and 3 body 
states with comparisons to various bound state equation results in x^0 theory 
is currently under study and will be presented in a separate article. |]19| 



In this paper we present a complete numerical algorithm for the evalua- 
tion of n-hodj masses {n < 3) and distribution probabilities for scalar x^0 
interaction in 3+1 dimensions. By providing this algorithm we intend to 
facilitate the comparison of various nonperturbative methods with the exact 
quenched results in x^0 theory. The organization of the paper is as follows: 
In the next section we present a brief summary of the FSR formalism. In 
the third section, using various 1, 2, and 3-body cases as examples, we dis- 
cuss how results are obtained. And in the fourth section we explain the 
components of the program. 



2 The Feynman-Schwinger representation for 
scalar fields 

We consider the theory of charged scalar particles x of rnass m interacting 
through the exchange of a neutral scalar particle of mass /i. The Euclidean 
Lagrangian for this theory is given by 

Ce = x* [m" -d^+gct>]x + \ 0(/^' - 5')0. (1) 

The two body Green's function for the transition from the initial state $j = 
X*(x)x(2;) to final state $/ = X*{y)x{y) is given by 

G{y,y\x,x) = N j Vx* j Vx j V(t)^)^,e~'^. (2) 

The final result for the two-body propagator involves a quantum mechanical 
path integral that sums up contributions coming from all possible trajectories 
of particles 

G = - P ds P ds I {Vz),y I {Vz)sye-^^^\ (3) 

where S[Z] is given by 

S[Z] = -iK[z, s] - iK[z, s] + Vo[z, s] + 2Vi2[z, z, s, s] + Vo[z, s]. (4) 

where 

2 , . , 1 /-i dz^ir) dz'^ir) 



K[z,s] = {m' + ie)s-— dr ^^ ^ / \ (5) 

^ ^ ^ ' As Jo dr dr ^ ' 

Vo[z,s\ = ^-^s' f'dr f'dT'A{z{r)-z{r'),f,), (6) 
/ Jo Jo 

(-,2 „i „i 

Vi2[z,z,s,s] = —ss dr df A{z{r) - z{f), /j,). (7) 
2 Jo Jo 

Here the Vo[2;,s] term represents the self energy contribution, while the 
V^i2[2;, ^, s, s] term represents the exchange interaction (Fig. [ip. The inter- 
action kernel A[x) is defined by 

A(x,/i) = / /^, f\ = -4r^KAi2\x\). (8) 

^ '^^ J (27r)>2 + ^2 47r2|a,| ^vpi \j \ ) 



While we present expressions for the 2-body case, generahzation to an arbi- 
trary n-body system is trivial. The bound state spectrum can be determined 
from the spectral decomposition of the two body Green's function 

oo 

G(T) = J2 Cne-™''^, (9) 

n=0 

where T is defined as the average time between the initial and final states 

T = -(y4 + l/4-a;4 -X4). (10) 

In the limit of large T, the ground state mass is given by 

While this result in principle is correct, the convergence to the asymptotic 
mass is slow due to the continuum contribution. The spectrum of the particle 
involves the mass pole and a cut beyond this pole representing the contin- 
uum contribution. Assuming that the constituents of the bound state are 
restricted to be at equal times in the initial and final states, the Green's 
function can be written as a function of time T, total displacement R, and 
relative coordinate r, G(T,R,r).0 

In order to eliminate the contribution of continuum states, introduce the 
Fourier transform, 

G{T, P,p)= f rf^R dh e'P-^+'P"G(T, R, r) (12) 

where P is the CM momentum, and p is the relative momentum between 
particles. Setting both P = p = one has 

G{T,0,0) = f d^Rdh G{T, R, r) , (13) 

which eliminates the contribution of the continuum and projects out the s- 
wave state. While an integration over r is not necessary for the elimination 
of the continuum contribution, it is useful in eliminating the contribution of 
states with nonzero orbital angular momentum. 

^Dependence on the initial relative final coordinate rp is implicit. 



While the result for the Green's function Eq. ^ is exact in the quenched 
approximation, due to its oscillatory behavior it is not appropriate for Monte- 
Carlo simulation. In Ref. [Q it was shown that one can perform a Wick 
rotation in variable s to avoid these oscillations. In the limit g^ ^ the 
dominant contribution to the integral in Eq. (^) can be shown, by using the 
saddle point method, to come from 

T 
s = iso = i-—. (14) 

2m 

Since the large s values do not contribute to the integral even without the 

interaction term, it is a good approximation to suppress the g^ term at large 

s values. While this suppression is done it is important that the integrand 

is not modified in the region of dominant contribution s ~ isq. This can be 

achieved by scaling the s variable, in the interaction term only, by 

s 



R{s,so) 



(15) 



where 



R{s,So) = l-{s-tSo)yT^. (16) 

In the free case, {g^ = 0), the width W of the region of dominant s contri- 
bution goes as 



Therefore, in the free case the dominant contribution to the s integral comes 
from i{so — W) < s < i{so + W). This claim is supported by the Monte-Carlo 
simulation results. In Fig. ^ we present the results for s-distributions for two 
different coupling strengths for time T = 40, and m = 1 GeV. According 
to the estimate given above Eq. (P^D, for g = 0, the dominant contribution 
to the s-distribution comes from the region 15.53 < s < 24.47 which is in 
agreement with the result presented in Fig. ^j. In order to ensure that the 
scaling given in Eq. (|1^) does not make a significant change in the region of 
dominant contribution, F should be chosen such that 

T>W (18) 

As one increases the coupling strength, the value of Sq deviates from its free 
value. Therefore, in general, Sq has to be defined self consistently by monitor- 
ing the peak of the s distribution. In Figure ^ we display the s-distribution 
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for two different coupling strengths. It is seen that as the couphng strength 
is increased the peak of the s-distribution moves towards higher s values. In 
general the peak of the distribution can be parameterized as 

The dependence of C on coupling strength g^ is determined self consistently. 
In Fig. H C, which gives the location of the stationary point through Eq. ([T9|), 
is plotted as a function of the coupling strength g"^. According to Fig. ^, it 
is not possible to find a self consistent stationary point beyond the critical 
coupling strength of g^ = 31 GeV^. A similar critical behavior was also 



observed in Refs. [Q within the context of a variational approach. 

The insensitivity of the dressed mass to the width F has been inves- 
tigated and found that a choice of F^ = 2W'^ was satisfactorily large. 
Results presented here employ the same value of F. 

Prescription given by Eqs. (p!5[), and (0) enables one to perform a Wick 
rotation in variable s to obtain a finite and non-oscillatory expression for the 
fully interacting two-body propagator: 

G = J dsj ds J iVz),y J {Vz)^y (20) 

K[Z, S] - K[Z, S] + Vo[z, Sr] + Vo[z, Sr] + 2^12 [-2, Z, S^, Sr 



xexp 

where 



(21) 



R{s,so) 
The path integral is discretized using 

{VU - {N/iTisf'ul-,' J d\, (22) 

where the s-dependence is critical in obtaining the correct normalization. 
Note that the integration over the final coordinates is not included in this 
expression. 

The discretized versions of kinetic and interaction terms are given by 

N ^ 
K[z,s] -^ {m^ + te)s-—J2{zi-z,.if, (23) 



g^s^ ^ 1 
Vo[z,s] ^ -—Y.A{-{zi + Zi_i-Zj-zj_i),fi), (24) 

Vi2[z,z,s,s] -^ ^^ ^ A(-(z, + 2;i_i-Zj--Zj-_i),/i). (25) 

We next address the regularization of the ultraviolet (short distance) 
singularities. The ultraviolet singularity in the kernel A(x,/i) Eq. (P) can 
be regularized using a Pauli-Villars regularization prescription. In order to 
do this one replaces the kernel 

A{x,fi) — > A{x,fi) - A{x,afi), (26) 

where a is in principle a large constant. The ultraviolet singularity in the 
interaction is of the type 

j dzzA{z,ii). (27) 

At short distances the kernel A{z,^) goes as l/z"^. Therefore, the self en- 
ergy calculation involves a logarithmic type singularity. The Pauli-Villars 
regularization takes care of this singularity. The Pauli-Villars regularization 
is particularly convenient for Monte-Carlo simulations since it only involves 
a modification of the kernel. In order to achieve an efficient convergence 
in numerical simulations we use a rather small cut-off parameter a = 3. 
This choice leads to a less singular kernel. The value of a can be increased 
arbitrarily at the cost additional computational time. It should be noted 
that while we employ a Pauli-Villars regularization throughout this paper, 
in general, the bound state problem without self energies does not have any 
UV singularities. UV singularities are only associated with the self energies of 
the particles. After this brief summary of the formalism, in the next section 
we present some of the results and discuss the details of the algorithm. 

3 Running the code: applications 

The Monte-Carlo simulation starts by choosing an initial configuration for the 
trajectories of the particles. The choice of the initial trajectory is arbitrary 
(except at the end points). In Fig. ^the evolution of the action as a function 
of number of updates for two different initial conditions is displayed. Starting 



with a random initial trajectory is analogous to a high temperature system 
and the termalization (reaching to the ground state) takes a long time (about 
1000 updates). However if one starts by a classical free trajectory, the initial 
configuration is analogous to a frozen system without any fluctuations and 
the termalization results in an increase of the action. However asymptotically 
both results should converge as shown in Fig. ^ We usually start with an 
orderly system and disregard the first 1000 updates. Step sizes of the random 
walker in configuration space and in s-space should be chosen such that the 
average acceptance rate of Monte-Carlo updates are about % 50. In each 
run we typically make 500000 updates of trajectories. In order to reduce the 
statistical errors runs must be repeated (usually more then 10 times) using 
different random number seeds. The correlation function X{n) of sampled 
configurations in each run is defined as 

{m) 

where m{i) is the mass measurement at the i'th update. The correlation 
function X{n) measures how the information about a given configuration 
is lost as a function of the number of updates n. In Figure ^ we show 
the correlation function for the 1-body problem. According to Fig. ^ the 
number of updates necessary for the correlation between sampled trajectories 
to vanish is around n = 1000. The number of configurations sampled in our 
simulations, which is around 500000, is well above the correlation length 
1000, thereby insuring that uncorrelated trajectories are sampled during the 
Monte-Carlo integration. 

The time T required to reach to the asymptotic limit increases as g"^ 
increases. For g"^ = 25 the asymptotic value of the mass, given by Eq. ([TT|) , 
is obtained around mT = 40. In particular for the result shown in Fig. ^, as 
g"^ ranged from g"^ = GeV^ to (7^ = 31 GeV^ the asymptotic time values 
used were increased from mT = 35 to mT = 45. As time T is increased 
the number of steps A^ should also be proportionally increased so that the 
step size of the particle trajectory remains the same. As one increases the 
coupling strength g"^, trajectories of the particles deviate from the classical 
trajectory to a greater degree. This increase in fluctuations requires that one 
uses a higher A^ value. The number of steps A^ particles take between the 
initial and final coordinates is typically chosen to be 35 < A^ < 45. 
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In Fig. I the g"^ dependence of the one-body dressed mass is presented. 
1-body masses presented here are lower than those in |^. This is due to the 
fact that in that work C was assumed to be approximately 1. However it 
was subsequently reahzed that C deviates from 1 significantly (see Fig. |^) 
as coupling strength is increased, and a self consistent determination of C 
is important to pin down the critical point. According to Fig. ^ there is no 
self consistent stationary point Eq. (0) beyond the critical coupling strength 
(?^ = 31 GeV^, and the 1-body dressed mass becomes unstable for (yf^ > 31 
GeV2 (Fig. I). 

Next two applications involve the two and three-body bound states of 
equal mass particles. While the algorithm provided is capable to take into 
account all self energy and vertex dressing corrections to the bound state, in 
the bound state applications provided here the self energy contributions of 
particles were not taken into account. This is controlled by a switch in the 
input file (see Table ^]1| for input options). In Fig. ^ the g"^ dependence of 



the two-body bound state mass is presented. Beyond the critical coupling 
strength of g^ = 100 GeV^ the 2-body mass becomes unstable. During 
the Monte-Carlo simulation the radial distributions of particles are stored in 
histograms. In general for an n-body bound state there are n(n— 1)/2 relative 
distances. Since the particles are assumed to have equal time coordinates in 
the final state, essentially the relative distance is equal to the spatial distance. 
Let rij = \ri — rj\ be the distance between particles i and j. In this case the 
histogram Pij{rij) stores the number of final state configurations sampled in 
the interval (r^j — 6/2,rij + S/2), where 6 is the bin size. The range of r^j 
values is specified in the program using 

20 , , 

< Tij < Fermi. (29) 

nii + rrij 

The choice of range is arbitrary as long as a reasonably smooth histogram 
is produced. The number of bins in each histogram is chosen to be 100. 
Therefore for two particles of mass 1 GeV, the maximum radius to be his- 
togrammed is 10 Fermi; and the size of each bin is given by 5 = 0.1 Fermi. 
In Fig. Ijwe present the projection of the radial probability distribution onto 
a two-dimensional surface. The two-body probability distribution, presented 
in Fig. ^, shows the result for g^ = 25 GeV^. Notice that the probability 
distribution vanishes at the origin. This is due to the phase space factor 
47rr^. In order to be able to compare the probability distribution with a 
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wave function the phase space factor should be factored out. In general the 
radial wave function amplitude \l/(r) is given by: 




^(0l = \7±i (30) 



The normalization of the probability distribution histogram presented in 
Fig. ^ is arbitrarily fixed such that the maximum entry is equal to 1. 

In generating the surface plot (Fig. |), it is assumed that one of the parti- 
cles is fixed at the origin. The amplitude of the surface gives the probability 
of finding the second particle at a distance r. While the surface plot for 
a two-body bound state is not necessary for a visual understanding of the 
bound state structure, the three-body bound state demands a three dimen- 
sional plot. In the three-body case there are three relative coordinates. In 
order to be able to plot the three-body probability distribution we fix the 
location of two of the particles along the y axis, and calculate the probability 
of finding the third particle in an arbitrary distance from two fixed particles. 



In Fig. |T0| we represent the probability distribution of the third particle for a 
given fixed configuration of the first and second particles. Assume that the 
fixed particles are particle 1 and particle 2. The probability distribution of 
the third particle is given by 

Psdrsl) = Pisilrs - riDPssdrs - r^\) (31) 

Probability distribution PsdrsI) includes the phase space contribution. There- 
fore PadrsI) represents the probability of finding the third particle in a ring 



shown in Fig. |Tl|. For example in the first plot of Fig. |T0| two fixed particles 
are very close to each other such that the third particle sees them as a point 
particle. Therefore the probability distribution of the first plot of Fig. [1^ is 
very similar to the two-body distribution given in Fig. ^ However as the 
fixed particles are separated from each other the third particle starts having 
a nonzero probability of being in between the two fixed particles (second and 
third plots of Fig. ^). Eventually when the two fixed particles are kept away 
from each other the third particle has a nonzero probability distribution only 



at the origin (the last plot shown in Fig. |T^). In Fig. 12 the two-body dis- 
tribution function P12 in a three-body system is shown. All of these plots 
were produced with three equal particles of mass 1 GeV, and the coupling 
strength of g"^ = 64 GeV^. 
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4 Program details 

In this section we discuss the details of the program. We start by providing 



the tables of input parameters ^^ and arrays [4.2| and then summarize the 
components of the program. 



4.1 Description of Input parameters 



SIG = 1 GeV 


Arbitrary momentum scale 


IPAR(I) 


1=1,2,3, Particles present (0=n, l=y) 


IEXCH(I,I) 


1=1,2,3, Self Interactions (0=n, l=y) 


IEXCH(I,J) 


(/, J) = (1,2), (1,3), (2, 3) Exch. Int. (0=n, l=y) 


QM(I) 


Mass of particle I 


G 


The coupling strength 


XMU 


The exchange mass fi 


ALPHA 


The Pauli-Villars mass ratio: rripjj = a fi 


BETA = 1 


The width of R{s, sq) = 1 - (3{s - iChf/T"< 


GAMMA = 1 


The T dependence of the width. 


C 


The peak location of the s distribution s=Csq. 
C is determined self-consistently. 


N 


Number of steps along the trajectory 


NSMPL 


# of sampled trajectories in MC integration 


NVOID 


No. of uncounted samples for initial termalization 


ZSTEP 


Max. step size of the random walker in coord, space 


SSTEP 


Max. step size of the random walker in s space 


EPSA 


The short distance(UV) cut-off for preventing 
exphcit occurrence of 1/0.0 type of singularity. 
This is not a regulator however. The UV 
regularization is done using Pauli-Villars subtraction 


IDUM 


The seed of the random number generator 


Z(I,0,J) 


Initial coordinates of particle 1=1,2,3 


Z(I,N,J) 


Final coordinates of particle 1=1,2,3 


IWRTA,IWRTM 


Action, mass to be stored in file ? (0=n, l=y) 


INTOUT 


Integrate over final coordinates ? 
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4.2 Arrays 



QM(3) 


Particle masses 


IPAR(3) 


Determines which particles are present. 


Z(3,0:400,4) 


The trajectories of particles 


ZNEW(3,0:400,4) 


The updated trajectories of particles 


SMAX(3) 


The maximum value of s value to be histogrammed, 
(not an integration cut-off) 


RMAX(3,3) 


The maximum relative distance between particles 
to be histogrammed, (not an integration cut-off) 


SHIST(3,200) 


Histogram of s values 


RHIST(3,3,200) 


Histogram of relative distances between particles 


SUMK(3) 


Kinetic energies of particles 


SUMV(3,3) 


Potential energies of particles: SUMV(I,I): self 
energy of the I'th particle, SUMV(I,J): the exchange 
energy between I'th and J'th particles. 


SUMKN(3) 


Updated kinetic energies 


SUMVN(3,3) 


Updated potential energies 



4.3 Main Program 

The main program starts by calling the INPUT subroutine. This subroutine 



reads input parameters described in Table ^^. Next, subroutine XINTGR 



is called to perform the Monte-Carlo simulation. In the following the role of 
each subroutine and function is explained. 

4.4 Subroutine INPUT 

In this subroutine 12 input parameters are read. In addition to reading 
these parameters, histograms for the radial and s distributions, and initial 
trajectories of the particles are initialized. Initialization of trajectories is 
done using the classical free trajectories of particles. 
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4.5 Subroutine XINTGR 

This subroutine controls basic steps of the Monte-Carlo samphng process. 
First, the kinetic and the potential sums corresponding to the initial tra- 
jectories of particles are calculated by calling the XSUMS subroutine. Next 
step is the termalization process. In order to reach termalization the con- 
figuration of trajectories are updated NVOID times by calling the UPDATE 
subroutine. First NVOID updates are not used in the actual calculation of 
the bound state mass or the probability distribution. After the termalization 
the sampling is done for NSMPL times. The results for the s values and the 
relative distances of particles are histogrammed. Finally the bound state 
mass is calculated. 

4.6 Subroutine UPDATE 

This subroutine is responsible for updating the current configuration of tra- 
jectories. Each coordinate and s parameter is updated once. Samplings are 
done according to distribution 

e-^W (32) 

distribution. If the ratio r, 

r = e-(^M-^W\ (33) 

is larger than 1 updates are always accepted. When r < 1 the update is 
accepted with a probability of r. 

4.7 Subroutine XSUMS 

This subroutine calculates kinetic and potential sums before the configuration 
updates start. Since calculation of the kinetic and potential sums is a costly 
operation, after every update we calculate the shift in the sums. 

4.8 Subroutine ACTION 

This subroutine calculates the action using known values of kinetic and po- 
tential sums. The action is stored in a file to monitor the termalization. ( 
See Fig. I). 
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4.9 Subroutine UPDSUM 

Since calculation of the kinetic and potential sums is a costly operation, after 
every update we calculate the shift in the sums. Using this shift sums are 
updated. 

4.10 Function XDERIV 

This subroutine calculates the derivative operator S'[Z] Eq. ( [Tl| ) which gives 
the mass of the bound state. 

4.11 Function DELTA, DELTAP 

These subroutine calculate, respectively, the interaction kernel Eq. (H) and 
its time derivative. 

4.12 Function DLARAN 

This is a random number generator obtained from LAPACK package at 
Netlib. |]T^ DLARAN returns a random real number from a uniform (0,1) 



distribution. In the actual implementation of the program we have used a 
Numerical Recipes random number generator (ran2). However for copyright 
reasons here we provide this alternative random number generator. 

4.13 Functions BESKO, BESKl, KZEONE 

BESKO and BESKl are modified bessel functions Ko{x) and Ki{x) which are 
needed in calculation of the interaction kernel and its derivative. BESKO 



and BESKl call KZEONE subroutine |T8|. KZEONE subroutine returns 
real and imaginary parts of e^Ko{x) and e^Ki{x). The KZEONE subroutine 
is considerably slower than the Numerical Recipes subroutines for modified 
bessel functions. For copyright reasons we provide KZEONE rather than the 
Numerical Recipes subroutines for the modified bessel functions. 
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Figure 1: A sample trajectory of each particle along with various interactions 
are shown. Except the matter loops all self energy, exchange energy and 
vertex dressings are included. 

5 Conclusions 

In this paper, using the Feynman-Schwinger representation, we have pre- 
sented an algorithm that calculates 1,2,3 body masses and distribution prob- 
abilities in the quenched approximation for x^0 theory in 3+1 dimension. 
The FSR approach provides an efficient method to calculate nonperturba- 
tive propagators. In this work we have presented results of applications to 
the 1, 2 and 3 body states. A detailed comparison of quenched bound state 
results of the FSR approach with the bound state equation predictions is 
under study and will be presented in a separate physics article |T^. It is 
hoped that, through comparison, this simple and rigorous nonperturbative 
method will enhance our understanding of various nonperturbative bound 
state models and approximations. 
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6 Test run INPUT 



1.0 


SIG (GeV) 


100 


IPAR(l), (2), (3) 


100 


IEXCH(1,1)(2,2)(3,3) 


00 


IEXCH(1,2)(1,3)(2,3) 


1.0 


QM(1) 


1.0 


QM(2) 


1.0 


QM(3) 


2.0 


G 


0.15 


XMU 


3.0 


ALPHA 


0.5 


BETA 


1.0 


GAMMA 


1.0 


C 


35 


N 


500000 


NSMPL 


5000 


NVOID 


2.1 


ZSTEP 


4.5 


SSTEP 


l.OD-4 


EPSA 


1 


IDUM 


0. 0. 0. 0. 


Initial coordinates 


1. 1. 0. 0. 


u 


-1. 1. 0. 0. 


u 


0. 0. 0. 35. 


Final coordinates 


1. 1. 0. 35. 


u 


-1. 1. 0. 35. 


u 


00 


IWRTA, IWRTM 


1 


INTOUT 
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7 OUTPUT 



SIG 


1.0 GeV 


Particles present 


1 00 


Self interactions 


1 00 


Exchange " 1-2 





" " 1-3 





" " 2-3 





QM(1) 


1. GeV 


QM(2) 


1. GeV 


QM(3) 


1. GeV 


G 


2. GeV 


XMU 


0.15 GeV 


ALPHA 


3. 


BETA 


0.5 


GAMMA 


1. 


C 


1. 


N 


35 


NSMPL 


500000 


NVOID 


5000 


ZSTEP 


2.1 1/GeV 


SSTEP 


4.5 l/GeV2 


EPSA 


O.lOE-03 GeV-i 


IDUM 


1 


Initial coordinates: 


GeV-i 




110 GeV-i 




-110 GeV-i 


Final coordinates: 


35 GeV-i 




1 1 35 GeV-i 




-1 1 35 GeV-i 


INTOUT 


1 


z update % 


50.33 


s update % 


53.63 


Bound state mass 


0.97±0.001 GeV 
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0.9 



X 



10 



20 
g' (GeV') 
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Figure 3: The dependence of the peak of the s-distribution on the couphng 
strength is shown. The peak location is given by Sq = CT /2m. Beyond the 
critical coupling strenght of g^ = 31GeV^ a self consistent determination of 
C is not possible. Therefore beyond the critical coupling strength 1-body 
mass becomes unstable. 
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2000 



1500 



I 1000 
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500 



random initial trajectory 

classical initial trajectory 






''•4AiMWIit#|teM|Mkai|u||^ 



10 100 1000 

number of updates 



10000 



Figure 4: The termalization for two different initial configurations is shown. 
Irrespective of the initial conditions the action approaches to an asymptotic 
limit. The first fOOO updates are ignored. 
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X 



0.8 



0.6 



0.4 



0.2 



-0.2 



Correlation function 



■^■^^»,.>^j"'^'Sv'VW"^i"VVy^V^ — 



10 100 

n (number of updates) 



1000 



Figure 5: Here we present the correlation function for the 1-body simulation. 
The number of updates necessary for the correlation to vanish is around 
n = 1000. The typical number of updates performed in our simulations, 
which is around 500000, is well above the correlation length, thereby insuring 
that statistically uncorrelated configurations are sampled. 
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1.1 



1.0 - 



0.9 



> 

CD „ „ 

O 0.8 



0.7 



0.6 



0.5 



0.0 



1-body mass 



I 



I 



I 



I 



10.0 



20.0 
g'(GeV') 



30.0 



40.0 



Figure 6: The coupling constant dependence of the 1-body dressed mass 
is shown. The asymptotic time value T used to obtain these mass values 
increase as g increases. While at g^ = GeV^ mT = 35 is sufficient, at 
g'^ = 31 GeV^ mT = 45 is needed. The peak of the s distribution was self 
consistently determined. Beyond the critical coupling strength of g^ = 31 
GeV^ the 1-body mass becomes unstable. 
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2.0 



2-body mass 



> 1-9 



1.8 



1.7 I • ' • ' • ' • ' • 1 

20.0 40.0 60.0 80.0 100.0 120.0 

g' (GeV') 



Figure 7: The coupling constant dependence of the 2-body bound state mass 
is shown. Beyond the critical coupling strength of g^ = 100 GeV^ the 2-body 
mass becomes unstable. 
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Figure 8: The 2-body probability distribution when the 1st particles is at the 
origin. The vanishing of probability at the origin is due to the phase space 
factor Eq. (BOl). 
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Figure 9: The radial probability distribution for the 2-body bound state is 
shown. 
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Fermi 



Figure 10: Evolution of the probability distribution for the 3rd particle is 
shown as the distance between the two fixed particles is increased. When 
the fixed particles are very close to each other the third particle sees them 
as a point particle (the upper left plot). As the fixed particles are separated 
from each other the third particle starts penetrating between them (2nd and 
3rd plots), and as the two fixed particles are maximally separated the third 
particle spends most of its time in between the two fixed particles (the lower 
right plot). 
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Figure 11: The radial probability distribution for the 3-body bound state 



shown in Fig. |iy includes the phase space factor represented by the ring in 
this figure. Distances ri3, and r23 are fixed on the ring and each point on the 



ring contribute to the probability distribution shown in Fig. |10 
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Figure 12: The radial two-body probability distribution P12 for the 3-body 
bound state is shown. This distribution is obtained by integrating out the 
location of the 3rd particle. 
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